Patterns in Flowing Sand: Understanding the Physics of Granular Flow 
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Dense granular flows are often unstable and form inhomogeneous structures. Although significant 
advances have been recently made in understanding simple fiows, instabilities of such flows are 
often not understood. We present experimental and numerical results that show the formation of 
longitudinal stripes that arise from instability of the uniform flowing state of granular media on a 
rough inclined plane. The form of the stripes depends critically on the mean density of the flow 
with a robust form of stripes at high density that consists of fast sliding plug-like regions (stripes) 
on top of highly agitated boiling material — a conflguration reminiscent of the Leidenfrost effect 
when a droplet of liquid lifted by its vapor is hovering above a hot surface. 



Granular flows are ubiquitous in the environment and 
in industry but there are still no known equations for 
general granular systems. For flow on an inclined plane, 
however, progress has been made [U El El H] , and current 
theories give a reasonable explanation of uniform flow- 
ing states. The simplicity of the inclined- plane geometry 
plays a crucial role in such descriptions because boundary 
effects, which can be exceedingly complicated in granu- 
lar flows, are well deflned. This system is also well suited 
to Discrete Element Method (DEM) simulations 17] 
because it can be treated periodically and steady flows 
often result, so good statistics can be obtained. The com- 
bination of detailed simulations and experiments has led 
to a solid understanding of steady flow states in this sys- 
tem where experiments and simulations can be well sum- 
marized 1^ |3] [SI in] by a recently proposed model [8]. 
Little is known, however, about the instabilities of this 
model or indeed of most granular flows, and the model 
has only been well tested in simple shear states. In con- 
trast, the Navier-Stokes equation of fluid flow has been 
known for over a century, and its accuracy has been re- 
peatedly tested by comparing the results of linear and 
weakly-nonlinear stability analysis to experimental sys- 
tems displaying an instability from a simple state to one 
with a distinct pattern [TDj. Such an approach, that is, 
investigating the growth of instabilities from their respec- 
tive steady states [TTJ [T2], will certainly be very useful 
in testing granular constitutive models and will provide 
critical tests for emerging theories of granular flow. 

The steady and fully developed state of a rapid, dilute 
granular flow on a rough inclined plane was shown exper- 
imentally to be unstable to the formation of longitudinal 
vortices observed as lateral stripes [11]. In this pattern 
the downstream velocity and the layer height vary peri- 
odically across the flow consisting of higher-slower and 
lower-faster regions. The development is attributed to a 
mechanism analogous to the Rayleigh-Benard instability 
in heated liquid layers [13] • The average packing fraction 
rjav of this flow was below 0.3 corresponding to a relative 
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density pr = rjav/rjs ~ 0-5, where ijg is the static pack- 
ing fraction (rjs is in between the random-closed-packed 
and random-loose-packed packing fractions [2]). This 
low density state, however, is hard to find either numer- 
ically or for some materials experimentally. Instead, we 
show that when increasing the plane inclination angle, 
the stripe state that emerges naturally is an instability 
of a dense uniform fiow state, that this stripe state is ro- 
bust and easy to find, and that the maxima of the down- 
stream surface velocity correspond to the highest regions 
of the modulated height profile in qualitative agreement 
with the fiow rule for the uniform state [TS] . 

The fiow was analyzed for a wide range of granular ma- 
terials including different sized sand and glass beads and 
various copper samples with different particle shape. We 
demonstrate, using the apparatus illustrated in Fig. [T] 
that longitudinal stripes are robustly observed over a 
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FIG. 1: Schematic illustration of the experimental setup. 

broad range of flow conditions with relative densities 
in the range 0.2 < pr < 0.95 (corresponding to about 
0.12 < rjav < 0.57), separated into the dense flow state 
for Pr ^ 0.6 (i.e. rjav ^ 0.36) and the dilute stripes for 
lower Pr as described below. 

In the present study two experimental setups were 
used. The flrst apparatus, described in detail else- 
where |16j . was used to characterize flow regimes over 
a wide range of 6. Because the system was enclosed in a 
cylindrical tube, precise measurements of the height pro- 
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file h{y) and determination of the flow structure were 
performed in a second setup. It consisted of a glass 
plate with dimensions 2.27 m x 0.4 m inclined at an angle 

9 = 41.3°, see Fig. [T] One layer of the d — 0.4 mm sand 
grains was glued onto the surface of the glass plate in 
both setups to provide a random roughness. Measure- 
ments were repeated using sandpaper with characteristic 
roughness of 0.19 mm (grit 80) and produced very simi- 
lar results. The surface velocity of the flow, with down- 
stream and transverse components u"^ and w**, and the 
height profile h(y), determined by a laser sheet, were ob- 
tained simultaneously using two cameras (Fig.[lJ. Using 
camera 3, the velocity at the bottom of the layer (u^ 
and v'') was taken at a location where the glass plate 
was clean. The flow velocities were determined using 
particle image velocimetry on image sequences taken at 
2000 frames/s. The relative density pr was measured us- 
ing a method described in detail elsewhere [16,. The 
majority of the data presented in this paper was ob- 
tained with sorted sand with d = 0.4 ± 0.05 mm and 
d = 0.2 ±0.05 mm. The stripe state was also detected for 
glass beads with d — 0.18±0.05 mm, d = 0.36±0.05 mm, 
and for four sets of copper particles with similar size 
(d = 0.16 ± 0.03 mm) but various shapes. The angle 
of repose 9r varied between 20.9° < Or < 33.8° for the 
materials tested. 

DEM simulations were also performed to investigate 
the instability. A soft particle model was used with a 
damped linear spring for the normal force (coefficient of 
restitution 0.8) and Coulomb friction for the tangential 
force (coefficient of friction 0.5). Particle stiffness was 
chosen so that the maximum overlap was less than 1%. 
The time step was 1/10 of the binary collision time. The 
base was made of identical particles held at fixed posi- 
tions taken from another simulation where a thick layer 
was allowed to form randomly. All quantities were non- 
dimensionalized using the particle diameters and gravity. 
The instability was found over a range of parameter val- 
ues: slope angle 34-39°, restitution 0.80-0.95 and width 
greater than 50. Below we present results from one typi- 
cal simulation with a slope angle 37° . If the slope angle in 
the simulations is reduced below 9r then the flows come 
to rest with a packing fraction of 0.6. We use this pack- 
ing fraction to normalize the results. The system was 
periodic in the x direction (down-slope length 24.3) and 
the y direction (cross-slope width 120.15). The number 
of particles simulated was 55 761, so the volume of the 
particles over the xy area corresponded to a height of 

10 (height at 100% packing fraction). The system was 
run for several months until steady state was achieved. 
Contact stresses were calculated according to the method 
of [T7| using a Dirac delta function as the weight func- 
tion. The stresses, particle positions, particle velocities 
(first and second moments), were gridded with spacing 
0.1 vertically and 0.45 laterally using linear interpolation. 
These were calculated every program step and averaged 
over 500 time units. 

In the experiments the stripe state for all the mate- 



rials tested has qualitatively similar characteristics, and 
there is no sharp transition between the states observed 
in the dense and dilute flow regimes. Nevertheless, the 
flow structure of the dense flow state is quite different 
from the dilute flow case reported earlier pjj . The struc- 
ture of the dense stripe state consists of relatively-narrow, 
dense, fast-moving regions that are also the highest. In 
the dilute regime the fast moving region corresponds to a 
height minimum, as schematically illustrated in Fig.|2]D,c. 
In both cases, however, the grains sink in the fastest mov- 
ing regions. The continuous transition between the two 
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FIG. 2: (a) The phase diagram of the system presenting 
the mean flow density pr as a function of normalized plane 
inclination tan 0/ tan 6'^. and the normalized hopper opening 
H = H/d based on data taken for sand, glass beads and var- 
ious copper samples. Various flow regimes are indicated, the 
bullet in the dense stripe domain corresponds to the simula- 
tion data presented. Illustration of the flow structure in (b) 
dense and (c) dilute regimes, gray levels indicate local density. 

regimes is characterized by increasing height of the slow 
moving region as the plane inclination is increased as 
illustrated in Fig.[3^-c or in movies taken for various ma- 
terials at [H]. The density decreases with increasing 9 
in a similar way for all materials as shown in Fig. [2^, 
where pr is plotted as a function of the normalized hop- 
per opening H — H/d and normalized plane inclination 
tan 9/ tan Or- Generally, stripes are only observed for 
tan 6*/ tan 6*^ > 1.25. Stripes with the structure typical 
for the dense regime are observed for 0.6 < pr < 0.95 
(corresponding to 0.36 < rjav < 0.57), whereas the di- 
lute regime, when it is observed, exists in the range for 
0.2 < Pr < 0.7 (i.e. 0.12 < r]av < 0.42). In the fol- 
lowing, we characterize the stripe structure in the dense 
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regime using experimental data shown in Fig.[3]obtained 
for sand with d = 0.4 mm and d = 0.2 mm and numerical 
results shown in Fig. |4] obtained from DEM simulations. 

The flow pattern 18] has a downstream surface ve- 
locity u'^iy) and a lateral surface velocity The 
downstream velocity has a relatively large modulation of 
(^max ~ '^av)/^av — 0-2, whereas the lateral velocity is 
very slow with maximal value w^ax ^ 0-04 where 
denotes the average downstream surface velocity. Cross 
sections of the velocity of the fully developed state in 
simulations are shown in Fig. |4^,b showing very similar 




FIG. 3: Images of the flow taken in reflected light (illumi- 
nation from the right) and normalized downstream surface 
velocity u" = u" /^/gd as a function of the normalized trans- 
verse coordinate y = y/d for sand with d = 0.2 mm and down- 
stream distance from the outlet a; = 1.55 m at plane inclina- 
tions (a) e = 42.6°; (b) 48.5° and (c) 52.2°, corresponding 
to trniO / tanOr = 1.56, 1.92 and 2.19, respectively. Data ob- 
tained aX X = 2.13m for sand with d = 0.4mm at 9 = 41.3°: 

(d) height proflles h{y) taken at various hopper openings H, 

(e) laser-line intensity (exposure time 4 ms), (f) transmitted 
light intensity, and (g) dimensionless wavelength X/hav of the 
pattern as a function of the normalized mean flow thickness 
hav = hav/d. To adjust hav the hopper opening H was varied. 

downstream velocity modulation, but somewhat smaller 
lateral velocities. In the experiments the lateral velocity 
measured at the surface v^{y) and at the bottom of the 
layer v^{y), is of opposite direction in accord with the flow 
structure obtained from the simulations (see Fig.jljj). 

The experimentally observed height variation in 
Fig.jsji agrees nicely with the simulation results (Fig. |4^- 
d) . The fast stripes correspond to a higher narrow maxi- 
mum of the h(y) curve but another set of less pronounced 
maxima is present between them. Thus, instead of a si- 
nusoidal h(y) profile observed for the dilute regime, a 
more complex h{y) is seen where a higher, global maxi- 
mum corresponds to the fast flowing region and a lower, 
local maximum corresponds to the slower region. The 
double peak is also seen in the transmitted light inten- 
sity in Fig. [Sf, but the wavelength of other important 



measures of the pattern, e.g., the downstream velocity 
or velocity in the yz plane (see Fig. [3^-c and Fig. |4^- 
b) do not change during the dilute-dense transition so 
the emergence of a nonsinusoidal height profile does not 
correspond to a full frequency-doubling. The pattern is 
observed only in a finite range of the flow thickness with 
the strongest amplitude at 10 < /i < 18. The wavelength 
of the pattern A is related to the mean flow thickness hav 
as 2.8 < \/hav < 4.5 as shown in Fig. [3^, so the cross sec- 
tion of a roll is elongated as compared to the traditional 
Rayleigh-Benard rolls with nearly circular cross section 

Numerical simulations enable us to visualize spatial 
variations of the relative density. Fig. [4j;, and of the in- 
ertial number / in Fig. |4}i. The inertial number, usually 
defined for incompressible flows 0, can be extended for 
the compressible case. We deflne / = dy^\D'\/ ^/p, where 
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FIG. 4: Results of the DEM simulations: (a) cross section 
of the downstream velocity u, (b) speed in the yz plane (s = 
+ yj'^) with streamlines, (c) relative density pr — 77/O.6 or 
packing fraction -q, (d) the inertial number I, (e) dependence 
of relative density on / and best quadratic fit, and (f) depen- 
dence of effective friction /i on the inertial number /. Solid 
line is best cubic flt. Dashed line is best fit to Pouliquen 
model (Eq.2 in [8]). 

D' is the deviatoric strain tensor {D' = D — Tr(D)/3), 
p the density, p the pressure, and we use the norm 
\A\ = ^y^^iAA^'^)/2 . We do not consider normal pres- 
sure differences and define p = — Tr((T)/3, where a is 
the stress tensor and a' ~ a + p. This is equivalent to 
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the Pouliquen definition when the flow is incompressible 
(Tt(D) = 0). The inertial number is proportional to the 
shear rate and to the ratio of the collisional stress to the 
total stress. The flow has the highest density in the fast 
moving region where / (and shear rate) is lowest. We 
identify this region as a "plug" sliding fast on top of a 
"boiling" region with very low relative density and high 
inertial number and shear rate, a configuration reminis- 
cent of the Leidenfrost effect IH] when a droplet of liquid 
lifted by its vapor is hovering above a hot surface. Ex- 
perimental data visualizing the level of fiuidization at the 
surface (Fig. Is?) and the profile of the transmitted light 
intensity (Fig. 3f) fully agree with this picture. 

The relative density (or packing fraction) determined 
from the simulation data decreases monotonically with 
increasing shear rate, see Fig. [4^, and shows an amazing 
collapse over a wide range of densities and values of the 
inertial number /. This suggests that, at least for fast 
chute flows, a simple equation of state giving the pres- 
sure as a function of shear rate, and density is possible. 
To test the rheology we calculate the effective friction co- 
efficient /i by /Lt = -Tr{a'D')/p\D'\. This is the /Lt that 

minimizes the residual error cr' + • Simpler calcu- 

lations of /i, e.g., = -(Jxzjozz OT fj,= -axz/p, produce 
poor results (no collapse would be seen in Fig.jiJ^), due to 
the complicated strain fleld. This definition is an exten- 
sion of the Pouliquen model [S] to include compressible 
fiows. Fig. [4^ shows the effective friction coefficient as 
a function of / and demonstrates a reasonably good col- 
lapse. The data does not fit well with Pouliquen rheology 
and is much more strongly curved and appears to even 
decrease for large / (Fig. |4|^). At low shear rates /i in- 
creases with increasing I, but at / = 0.7 there appears 
to be a turnover above which fj, decreases with increasing 
I. We believe that this behavior of the system is a key 
feature leading to the instability. Namely, by increas- 
ing the fiow thickness above a certain value, the inertial 
number near the plane reaches a threshold above which 
the effective friction starts decreasing. As a consequence 
the local inertial number (shear rate) grows even more, 
leading to stronger fiuidization near the plane. At the 



same time the fiuidization drops in the upper layer (a 
plug develops) and this plug slides even faster on top of 
the expanded fiuidized region. This mechanism agrees 
with other results for chute fiows [5] , but in that case the 
simulations were too narrow for the lateral instability to 
develop. On the surface the plug absorbs material from 
the two sides as the surface fiuidization is larger in that 
region, so the plug grows. In the simulations, as the insta- 
bility develops the mean fiow speed decreases until a new 
lower velocity is reached at which point the instability has 
saturated. Thus, the instability may play an important 
role on steeper slopes, where no simple, steady state is 
expected, by increasing the effective viscosity. Though 
we have framed our discussion in terms of fi, it is diffi- 
cult to draw too many conclusions from these results for 
a Pouliquen rheology. Our simulation data shows density 
differences, normal stress difference and that the devia- 
toric stress tensor is not aligned with the strain tensor. 
These all disagree with the assumptions of the Pouliquen 
rheology indicating that a considerably more complicated 
rheology is necessary along with an equation of state to 
describe these flows. 

This system provides a very interesting case for study- 
ing granular rheologies because there is a complicated 
strain and stress fleld but very simple boundary condi- 
tions. Since the flow is steady, accurate measurements of 
all the flow variables in a simulation are possible, the only 
constraint being computer time. This system is therefore 
ideal for developing and validating granular theories in 
new ways. An intriguing possibility is that the lateral 
ridges and furrows observed in large rock avalanches [20] 
may be the result of the same instability. 
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